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ABSTRACT 


This report presents a formulation for the development of a finite 
element program for the elastic analysis of two-dimensional bodies using 
the eight-node isoparametric quadrilateral. The program solves for both 
plane stress and plane strain problems. 

A general development of the finite element formulation based on 
the isoparametric displacement functions is presented. 

The program structure is given in the form of flow diagrams with 
descriptions of the numerical procedure used to obtain the element 
stiffness matrix, and the solution method employed to solve for nodal 
displacements. 

Three numerical examples, a plate under uniaxial tension, a plate 
under pure shear, and a beam under pure bending are presented to illus- 
trate the capability and limitations of the element implementation. The 
first problem is solved exactly by the element, as predicted by the form 
of its displacement functions. However, in the other two problems the 
accuracy of the solution is highly dependent upon the slenderness of the 
element, the number of elements in the map, and the numerical integration 
•scheme used to build the element stiffness matrix. 

The report ends with some recommendations concerning map generation, 
simplification of the input data, and extensions to solution of plasticity 
problems. 


1 . INTRODUCTION 


The advantages of the Isoparametric elements over the well-known 
constant strain triangle (CST) have been substantially demonstrated in 
the last few years. Based on the experience we have gained from the 
CST, development of a finite element program for isoparametric elements 
is highly desirable. As is described below, among the advantages of the 
isoparametric elements is their ability to be coupled with other iso- 
parametric elements which may be convenient in many cases, and furthermore 
with other type of elements such as crack-tip elements currently being 
developed . 

Among the numberless members of the isoparametric family, the four, 
eight, nine, and twelve node isoparametric quadrilaterals are the most 
popular for two-dimensional analysis, not only because the high order of 
their displacement functions but for their relative low cost in terms of 
computer time and input preparation. The efficiency of any particular 
element type used will depend on how well the shape functions are capable 
of representing the true displacement field. At the time we initiated 
this work, very little was known about the differences in accuracy between 
the eight and nine node parabolic elements. Primarily because of storage 
limitations, in particular in the transition to the elastoplastic version 
of the program (currently in process), we chose the eight node quadri- 
lateral as the pattern element of the program. 

The four node flwadHTtfteral (figure 1.1a) provides a linear displace- 
ment distribution. It has proved to be an efficient element since the 
construction of the stiffness matrix can be done in closed form, while 



the higher order elements require the use of a numerical integration 
procedure. The nine and twelve node quadrilaterals {Figure 1.1c and 
l.ld) produce parabolic and cubic displacement distributions within 
the element, respectively. The latter, although is very accurate in 
the handling of a large variety of problems, presents two major draw- 
backs. First, it requires a 4x4 Gauss Integration procedure to build 
the stiffness matrix, which means both large storage and a high computer 
bill. Second, the bandwidth of the master stiffness matrix is usually 
larger than that of an eight node isoparametric for a comparable number 
of elements, which may provide storage problems, in addition to longer 
running time in the solving routine. 

These reasons in addition to the results presented here and 
elsewhere [1-4] prove the wise choice with the eight-node isoparametric. 

•k 

Numbers in brackets denote entries in the list of references following 
the text. 




2. OVERVIEW 


The basic step of any direct stiffness finite element analysis 
is the unique description of the unknown functions {6}, in this case they 
correspond to the displacement field within each element, in terms of n 
parameters {6 q } given by 

{6} = CN] {6 q } (1) 

where n = number of nodes for element, 

{6 q } = vector of nodal displacements, and 
[N] - matrix of shape functions evaluated at the point in the 
element where displacements are desired. 

With displacements known at all points within the element, the strains 
at any point can be determined, resulting in a relationship of the form, 

(c) * [B] {6} (2) 

where the matrix [B] is conformed by the partial derivatives of the shape 
functions given explicitly in Appendix A. 

To satisfy equilibrium it is required that the forces acting on an 
element lumped as {f} concentrated at the nodes must balance the stresses 
within the element. That is, 

/ [B] T {o}dA - {f} = 0 (3) 

A 

where A is the cross sectional area of the element where tractions are 
prescribed, and (a) is the stress vector. 

For linear elasticity, the constitutive law between stresses and 
strains is • 


where [D] is the elasticity matrix that differs slightly between the 
cases of plane stress and plane strain, and {a Q } and {e Q } are the 


initial stresses and strains, respectively, 
into (3) we obtain for any element that 

Substituting (2) and (4) 

tf) = [k]{6T 

(5) 

with 


[k] - / [B] T [D][B]dA 
A 

(6) 

Once equilibrium is established at each 

node, and nodal displacement 

continuity is ensured, we may write for the complete structure a load- 

displacement relationship of the form, 


{F} = [K] {6} 

(7) 


where [K] is an £x£ matrix, l being the total -number of degrees of 
freedom of the structure and is called the overall stiffness matrix. 


To select an isoparametric element is to select a determined set 
of shape functions, which is not arbitrary. These are three minimum 
conditions which must be satisfied [1] in order to ensure convergence 
of the solution to the correct results: 

(a) The shape functions must be continuous between 
elements; 

(b) In the limit, as the element size is reduced to 
infinitesimal dimensions, the shape functions must 
be able to reproduce a constant strain condition 
throughout the element. This means that the 
unknown function must be able to take in the 
limit any linear form throughout the element. 
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(c) Rigid body motion is accomplished if every shape 
function satisfies the relation, 


n 

z NU,,n,) = 1 
i=l 1 1 


2.1 Shape Functions for the 8-Node Quadrilateral 

In the isoparametric formulation, the general relationship between 
the global cartesian coordinates (x,y) and the local curvilinear 
coordinates (?,n) is (Figure 2.1): 


x = N^x-j + + ... = {N}^{x n > 

( 8 ) 

y = N,y, + N 2 y 2 + . . . = {NJ T ty n > 

where = f(5,n) are the isoparametric shape or displacement functions, 
and {x n }, {y n ^ are ^he column vectors of the cartesian coordinates. 

For any values of z and n the x and y coordinates can be found once the 
functions N are known. 

In finite element analysis of two-dimensional stress problems, it 
is necessary to define the variation of nodal displacement components 
u. and v.j in terms of the nodal values of the displacement functions N^. 
Thus we have: 


u (l»n) s N-j + NgUg + ... + N n u n = {N}"^{u r| } 

v(S,n) = N-j v ] + N 2 v 2 + ... N R v n = {N} T {v n ) 

For the 8-node quadrilateral, suitable polynomials that describe the 
appropriate variation of the sides can be written as: 
x = a-j + :t 2 s + a 3 n + + a 6 n 2 


(9) 


2 2 
+ a-jZ n + ctga l 


( 10 ) 



or 

2 2 2 2 

x = [1 > £ » n> £ » Sn» n > 5 n* h s]{a^} 

which ensures that on the sides where n = + 1 (see Figure 2.1), the 

variation may be up to quadratic in s, and likewise when 5 = + 1 the 

variation may be up to quadratic in n* 

From (10) we can write: 

{X n ) = [C]{a n } 01) 

or 

■{«„> ■ [cr’ty 02) 

Or, in terms of the displacement functions we have: 

[N-j • N 2 , N 3 ...N b ] = [1, 5, n, 5 2 , n5, n 2 , ? 2 n, n 2 e][C ] _1 (13) 

Let ^ and be thecoordinate values of the i-th node in the normalized 
curvilinear system (s,n)> then the displacement functions for the 8-node 
isoparametric quadrilateral are given by: 


(a) 

For corner nodes: i = 1, 3, 5, 7 

* 


N. = 1/4(1 + 5^)0 + nn i )(?5 i + nn i - D 

04) 

(b) 

For midside nodes, 5. =0: i = 6 



N. = 1/2(1 - s 2 )(l + n^) 

(15) 

(c) 

For midside nodes, =0: i = 4,8 



N i = 1/2(1 + - n 2 ) 

06) 


2.2 Stiffness Matrix Formulation 

To assemble the master stiffness matrix of a structural system, 
the element stiffness matrices are properly superimposed. Thus, we can 
write the master stiffness matrix as: 

m 

[K] - e [kj 
i=l 1 


07 ) 


8 


where m = number of isoparametric elements and [k^] is the ith-element 

stiffness matrix defined by: 

[k,] = / [B] T [D]tB]dA (IB) 

1 A 

with A here defined as the element area. 

In the cartesian (x,y) system with corresponding displacement com- 
ponents (u,v), the components of strain in terms of displacements for 
linear elastic plane problems can be written as 


_ 3U 


and 


In matrix form we may write: 


‘ 3X 


. 3V 


' ay 


. 3U 

+ 21 

' ay 

+ 3X 

■3 

fi * 

3X 


0 

3 

ay 

3 

3 

ay 

3X J 


xy , 

Hence, for the element strain vector we write: 
formulation 


(19-a) 

(19-b) 

(19-c) 


( 20 ) 


for the isoparametric 


U) = £B]{u> 


where: 


[B] = C [B-j ] , [B 2 ],...oy] 


being m the number of nodal points of the element 
and.. 


{U> = tu, V-j U 2 V 2 





( 21 ) 


( 22 ) 


* » * 
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In terms of the displacement functions , the matrices can be 


written as: 



aN.. 


o 

FT 


o 

aN i 


W~ 

3N. 

8N i 

w 

FT 


( 23 ) 


Since N. = N.j(?,n), it is necessary to perform a coordinate transformation 
such that: 

3Nj 3N. 

Fx~ = T x(?,n) ’ slT ^ 

aN,. aN. 

' — T , 

sy “ y(s,n) sn 

where T, \ and \ are vector transformation functions relating 

yum/ 

x and 5, and y and n» respectively. 

In terms of the local coordinates we may write then. 
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rsx 


( 3N i ) 


|35 / 

3£ 


13X~f 

3X 

/ } S 



= CJ] 

i 

U 1 

3X 

ay. 

k l 

3N^ 

larj 

L 3n 

a n J 

lay 1 1 

ay 


(25) 


where [J] Is the Jacobian matrix, which can also be written in terms 
of the global cartesian coordinates as: 


U] « 


3N 1 

3N 2 


W~ ’ 

3T“”‘ 

* *3g 

3N-j 

3N ] 

3N 

in 

JrT ’ 

?n“ * * 

* *3n . 


x i ^ 

x 2 h 

i 


x i y 


(26) 


y y 

m 


where (x^>, {y^} are the global coordinates of the element nodes. 
Now we may write: 


/ aN i 


3N- \ 

3X 

■ mur 1 

iT 

3N., 

1 

3Ni ( 

3y 


3T/ 


(27) 


where I is the identity matrix, 


[I] 


.r °i 

Lo 1-1 


In equation (20) all the derivatives, 


3 N, 3N. 

5T and as 


can be obtained by hand from the original displacement functions N. 

for the particular type of element being used (see Appendix A). 

Note that independently of the number of element nodes [J] is always a 
2x2 matrix, for two-dimensional problems. 



Once (26) and (27) are calculated we may substitute after some 
manipulation into (25) and obtain the B^'s. Thereafter, the matrix [B] 
car, be assembled. 

Having [B] in terms of the normalized local curvilinear coordinates 
(5»n)> we must express the element of area in (18) in terms of the local 
coordinates % and n- 

Consequently, we write 

dA = dxdy = |J| d n d<- (28) 

and because we have used normalized local coordinates, the limits of 
integration are now -1 and 1 for both integrals. Thus, we obtain the 
element stiffness matrix in t^e form: 

[k,] » / / [BU.n)] T [D][BC 5>n )] |J| dfjdri (29) 

1 -I 

2.3 Numerical Integration Procedure . 

The difficulties of performing a. closed form area integration for • 
the element stiffness matrix (29) are avoided by use of 'numerical 
integration techniques. For elastic problems the closed /form integra- 
tion is tedious; for elasto-plastic problems, impossible, as in this 
case the constitutive matrix is not available in closed form. Before 
attempting any numerical integration, it is a requirement to 'know the 
degree of accuracy needed to ensure stability and convergence to the 
correct result. 

If the inti. element forces due to internal element stresses can be 
determined exactly by numerical integration for a constant state of 
stresses within each element, we can guarantee convergence to the- correct 
solution. This argument has been discussed in detail by Irons [6], 
Zienkiewicz [1], and others. 


The Interelement forces can be expressed as: 

tf) = (Jf EB] T [D][B]dA){u} (30) 

Or in terms of the stresses: 

Cf> = / [B] T {a}dA (31) 

A 

Thus, by numerical integration, we must be able to determine exactly 
such intergrals as 

jfJ /.! f i j i «*• 

obtained by combining equations (27) and (28). 

The integration using Gauss-Legendre quadrature is highly simplified 
because of the constant boundary values (-1, +1) that the local normalized 
curvilinear coordinates take in each element. 

According to Gauss-Legendre quadrature [5] we may write the element 
stiffness matrix given by (29) in the form: 
mm T 

[k] = E E C [ B ( ? ) !] [D][B(?^ ,0^)] 1 0 J (32) 

k 1 “1 

where are the coefficients of the integration points k and Z. 

T 

[B^k^)] and [BU^.n^)] are the matrices evaluated at the Gauss 
point (C^.n^), and m is the number of integration points used. 

In [1] it is argued that 2x2 Gauss quadrature provides an exact 
numerical evaluation of the element stiffness of the 8-node quadrila- 
teral, independent of the activated terms in the displacement functions. 
This is true in tension and shear problems where only terms as high 
as E? are present in (32) once the internal product is performed. It 
has been documented [7] that 2x2 Gauss quadrature is exact only when terms 


as high as cubic are present in the product form of (32). Several 
exercises performed during the development of ISOFINEL and its sub- 
sequent testing, showed that in problems such as pure shear, pure 
bending, three point bending, etc., with elements having a large 
slenderness ratio, much more accurate answers are obtained when using 
3x3 Gauss quadrature. 

The program in its actual version contains both features, at the 
user's option, both to facilitate the work discussed in this report and 
for possible subsequent usage. 


3. SOLUTION PROCEDURE 


The master stiffness matrix [K] in (18) is a square, symmetric and 
positive-definite matrix. The solution of the linear equations related 
to the finite element analysis takes a large fraction of the total CPU 
time. Furthermore, the storage space occupied by the complete matrix is 
sometimes so large that medium size problems are Insoluble if 
access to tapes or peripherial storage is not provided. However, at 
this point, we are interested in saving storage at the expense of time. 

The master stiffness in structural problems is banded when precau- 
tions have been taken in the nodal numbering scheme. Therefore, the 
storage of the complete matrix is unnecessary and is avoided by assembling 
K in a rectangular form at the same time it is constructed. The number 
of rows may correspond to the number of equations to be solved, and the 
number of columns equal to the maximum semi-bandwidth. 

Before procedi ng to solve the system of equations, the boundary 
conditions are applied. In the actual version, admissible boundary 
conditions are the following: 

- prescribed non-zero nodal displacements in the x-direction 

- prescribed non-zero nodal displacements in the y-di recti on 

- prescribed nodal forces in the x-direction 

- prescribed nodal forces in the y-direction 

- zero nodal displacements in the x or y direction. 

The program makes use of the Gauss elimination procedure [6] to 
solve for the nodal displacements of the expression 

(F> = [K]{6} 


4. PROGRAM STRUCTURE 


The ISOFINEL computer code has been developed having in mind two 
major premises. First, the main body of the program should remain intact 
once plasticity is implemented. For this reason, some parts have been 
partitioned into different subroutines that at first glance would look 
superfluous. Second, eventually the same program will host special 
crack-tip elements for the elasto-plastic analysis of stress fields at 
the vicinity of cracks. Thus, many parameters, such as number of element 
nodes, number of Gaussian integration points, etc,, have been left in 
terms of variables subject to be changed via input. 

In the following pages a flow diagram of the main program and its 
subroutines is presented. It is not the purpose here to expand into the 
particular aspects of the coding; some diagrams of subroutines whose 
contents are self-explanatory are omitted. 

4.1 ISOFINEL: Main Program 

4.1.1 Description of Terms : 

NEL = number of elements 

NPROB = number of problems 

ICONT = actual problem number (from 1 to NPROB), 

NRD = total number of degrees of freedom. 

4.1.2 Description of Subroutines Not Accompanied by Flow Diagrams: 

START: Reads material properties- geometry, dimensions, boundary 

conditions, etc., and print them. 

ZEROS: Initialize arrays for 

(a) Stiffness Matrix, [K] 


(b) Derivatives Matrix, [8] 

(c) Displacement Vector, {6} 

(d) Force Vector, {F > 

D I SPUN: Contains in explicit form the shape functions and their 
derivatives. Calculates them at Gaussian integration 
points for either the 2x2 G.I. option or the 3x3 G.I. 
STRSTR: Calculates the elements of the elasticity matrix, D, 
for either the case of plane stress or plane strain. 
BOUCO: Applies boundary conditions and modifies master stiffness 
matrix such that is in the correct form for the solution 
procedure. Possible boundary conditions are: 

- prescribed x or y nodal displacements 

- prescribed x or y nodal forces 

- zero nodal displacements. 

GAUSEL: Solves for the nodal displacements by using the Gauss 
Elimination Procedure. 

NODFOR: Performs the matrix multiplication 

[K] (6) 

to return the generalized nodal forces. 

OUTPUT: Prints out all information: 

(a) Nodal Displacements, (6) 

(b) Nodal Forces * {F} 

(c) Strains at Gauss Points, {e} 

(d) Stresses at Gauss Points, {a} 

(e) Stresses and Strains at the Elements Centroid 

(f) Principal Stresses 

(g) Coordinates (x,y) of Gauss Points 

(h) Execution time of Each Subroutine. 
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4.3 STIFEL: Stiffness Matrix Subroutine 
4.3.1 Description of Terms : 

NGA2 = 9 for 3x3 Gauss integration; 4 for 2x2 Gauss integration 
NEL8 = NEL*8 

NOEL = number of nodes for element (8, for present case) 

NDF = number of degrees of freedom per node (2, for 
present case) 

NBW = Bandwidth of master stiffness matrix 
INOEL = NRD*N0EL 

NM(I) = array of nodal configuration, (1=1, NEL8) 

XYM(J) = array of nodal coordinates, (J = 1, NRD) 

XJAC(I,J) = Jacobian matrix, (I = 1,2; J = 1,2) 

DDF (L ,K) = array of derivatives of displacement functions, 

(L = 1, NGA2; K = 1, INOEL) 

BM(I,J,K,L) = array for B-matrix, (1=1, NHL ; J = 1,2: 

K = 1, NOEL; L = 1, NGA2) 

TEMPK(I.J) - array for element stiffness matrix, (1=1, INOEL; 

J = 1, INOEL) 

STIF(K,L) = array for master stiffness matrix, (K = 1 , NRD; 

L = 1 , NBW) 



T 

MATMUL: Performs the matrix multiplication [B] [D][B] and returns 

the upper symmetric part of the element stiffness matrix. 



4.4 Flow Diagram 


ZQ 



ZERO ELEMENT 
STIFFNESS MATRIX 


ZERO JACOBIAN 
MATRIX, J 


CALCULATES THE 
JACOBIAN MATRIX, J 


© 






DO 20 JK = 1 ,NOEL 
KDI = NOEL* (KEY = 1) + JK 
K = NDF*(NM(KDI).- 1) + 1 
EPS(KEY,1,L) = EPS (KEY , 1 ,L ) 

+ BM (KEY , 1 ,JK,L)*DF(K) 
EPS(KEY,2,L) * EPS(KEY,2,L) 

+ BM(KEY,2,JK,L)*DF(K + 1) 
EPS(KEY,3,L) = EPS(KEY,3,L) 

+ BM(KEY,2,JK,L)*DF(K) 

+ BM(KEY,1 ,JK,L)*DF(K + 1) 

20 CONTINUE 


CALL STRESS (L) 


15 CONTINUE 


CALCULATES STRESSES AND 
STRAINS AT ELEMENT CENTROIDS 


10 CONTINUE 
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4.5 STRAIN: Strains and Stresses Subroutine 
4.5.1 Description of Terms : 

EPS ( I , J ,K) = Array of strains, e x ,e y ,y xy , (I = 1 ,NEL; 

J = 1.3; K = 1 ,NGA2) 

SIG(I,0,K) = Array of stresses, a x »V T *y* (I = 1,NEL; 
J = 1,3; K 1 , NGA2 ) 

DF(L) = Array of nodal displacements, (L = 1,NRD) 


4.6 Flow Diagram 



ZERO ARRAYS 
FOR STRAINS 
AND STRESSES 




:<x 


DETJ = XJAC ( 1 ,1 )*XJAC(2,2) 
- XJAC(1 ,2)*XJAC(2,1 ) 


DETERMINANT OF 
J MATRIX 


CHANGE «* XJAC (1 ,1 ) 

XU = XJAC (2 ,2 )/DETJ 
XI 2 « -XJAC (1 ,2)/DETJ 
X21 = -XJAC(2 JJ/DETJ 
X22 = CHANGE/DETJ 


CALCULATES THE 
ELEMENTS OF [J] 


DO 5 IJ = 1 , NOEL 
J = IJ*2 - 1 
K = J + 1 

BM(KEY,1 ,IJ,L) = XU*DDF(L,J) 
+ X12*DDF(L,K) 
8M(KEY,2, IJ ,L) = X21*DDF(L,J) 
+ X22*PDF(L,K) 

5 CONTINUE 


CREATES THE 
B-MATRIX 








5, NUMERICAL. RESULTS 


To test the capability and limitations of the program, two 
problems are solved: 

- Uniaxial Tension 

- Pure Shear 

- Pure Bending 

They are solved using both 2x2 and 3x3 Gauss integration. As it was 
described earlier, the displacement and strain solutions of these 
problems are contained in the shape functions of the eight-node iso- 
parametric quadrilateral. Hence, it is expected that they may be solved 
exactly by a single element. 

Three different maps containing 1,4 and 10 elements are tested. Each 
map is tested for three different element sizes. Figure 5.1 shows the 
configuration of each map. 

5.1 Uniaxial Tension 

The uniaxial tension problem is solved by prescribing uniform 
(/-displacements at the nodes on the surface y = L, and by preventing 
vertical motion of the nodes on the surface y = 0. Being the solution 
of this problem contained in the displacement functions of the element, it 
is not strictly necessary to prevent rigid body motion. However, it is 
done by restraining the middle node at the origin to move. For all maps, 
the prescribed displacements correspond to 1% of the original length of 
the plate. The solutions obtained for both displacements and stresses 
are in excellent agreement (less than 0.001% error) with the expected 
theoretical values. No differences are observed between the two integra- 
tion procedures studied. 







PURE SHEAR 

/r vs. Element Ratio, a/b 



Figure 5.2 




PURE SHEAR 

/x vs. Element Ratio, a/b 



Figure 5.3 






5.2 Pure Shear 


The pure shear problem Is solved by prescribing linear varying 
x-displacements at the boundary nodes. The x-displacements are given 
by the relation, 

«x = 

where is the displacement prescribed on the surface ij = L. For each 
map, <5 Q is 1% of the plate length, L. 

Several results are studied, Figures 5.2 and 5.3 show the variation 
of the normalized shear stress ,{r X y/T 0 ), with t Q being the theoretical 
stress, versus the element ratio (a/b),for the 2x2 and 3x3 Gauss integration 
procedures, respectively. The value of x chosen for each case is that 
at the Gauss point that differs more from the theoretical value. Values 
at Gauss points very close to the surface, where horizontal displacements 
are prescribed, are not taken in consideration because of local effects. 

These effects are only detected in Map 3, where Gauss points are closer to 
the boundaries, but they die out rapidly. Maps 2 and 3 present excellent 
agreement in both cases, being negligible the difference between the two 
integration procedures. However, for Map 3, Case C, there is a large 
•improvement when 3x3 G.I. is used, since for an element ratio of a/b = 1.0 
there is a 4.25% error in stresses in the 2x2 G.I. compared with a 0.001% 
error in the 3x3 G.I. When a/b = 2.0, the errors are 5 and 4.25%, 
respectively. Figures 5.4 to 5.7 present the variation of the normalized 
shear stress along the direction perpendicular to the plane of shear for 
a constant x/b, The value of x/b is determined by the position of the 
integration points within each element. In this case, x/b attains the 
values of 0.2887 and 0.3873 for the 2x2 and 3x3 G.I. procedures, respectively. 


PURE SHEAR 



Figure 5.4 




Figure 





PURE SHEAR 
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Figure 5.6 


PURE SHEAR 



Figure 5.7 







Maps 1, 2 and Case 3A of Map 3 provide solutions that are in excellent 
agreement with the existing solution. Neglecting the error introduced 
by the local effects at the plane of shear, the error in the shear stress 
for cases 3B and 3C is summarized in Table 5.1 below. 


Case 

CM 

X 

CM 

G.I. 

3x3 

G.I. 


Max. 

Error 

Mean 

Error 

Man. 

Error 

Mean 

Error 

3B 

1.4% 

0.2% 

1.5% 

0.2% 

3C 

2.0% 

0.6% 

3.0% 

0.6% 


Table 5.1 


It is observed that both procedures produce highly comparable results 
with a slight advantage of the 2x2 G.I. In this particular problem, 
the difference is due more to the numerical error introduced by the 
larger number of operations that take place in the assembling of the 
stiffness matrix in the 3x3 G.I. case. 

As it is known from the theoretical solution, the stresses o x and 
a are both zero for the case of pure shear. However in the numerical 
computation, the true value of these stresses is related among others to the 
element ratio, numerical integration scheme, and of course type of dis- 
placement functions used. In the present case we are interested only in 
the former two. Figures 5.8 and 5.9 are plots of the variation of the 
"amplified zero-stress" along the axis perpendicular to the plane of shear. 

The numerical values of a are chosen arbitrarily since the values obtained 

3 

for Oy are comparable. An amplification of 10 is chosen, since from the 
engineering viewpoint any pair of stresses with a difference in order of 
magnitude below this quantity is usually considered. From these two graphs 
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It Is appreciated the advantage of using 3x3 6.1. in more complex problems. 
Figure 5.8 shows for the 2x2 6.1, scheme, Cases 3A and 3B, that half 
distance away from the plane of shear the "assumed-zero stresses" start 
increasing in magnitude even that before being one-fourth of plate length 
close to the plane of shear have already attained values comparable to 
the shear stress. For Case 3C the same phenomena is observed well befrre a 
section half plate length away from the plane of shear. When the 3x3 6.1. 
scheme is used, this effect disappears in the first half section of the 
plate and is highly reduced in the rest. 

5.3 Pure Bending 

The pure bending problem is solved by applying linear varying tj - 
displacements on the surface y - L with a line of symmetry at x = 0. 

The prescribed displacements are given by 

*y = (*/b)5 0 

where <5 Q is the displacement prescribed at the corner nodes and corresponds 
to 1% of the plate length. Figures 5.10 and 5.11 show a variation of the 
normalized bending stress (o y /a Q ) with the element ratio for the 2x2 and 
3x3 6.1. schemes respectively. It is observed that by performing a 3x3 6.1. , 
the stress results improve enormously for Maps 1 and 2, but little is 
gained on Map 3. Figures 5.12 and 5.13 present the variation of the 
bending stress along the axial direction of the plate for Map 1. Although 
the maximum error in stress is only about 1.3% for the Case 1C of 2x2 6.1. , 
the error is reduced to 0.001% when 3x3 6.1. is used (Figure 5.13). Similar 
results are observed in Figures 5.14 to 5.17 for Maps 2 and 3. It is 
observed, however, that by moving to a 3x3 G.I., the improvement in the 
stresses is only of the order of 40%. These results are summarized in 
Table 5,2. 
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Case 

2x2 

G.I. 

3x3 

G.I. 


Max. 
Error % 

Mean 
Error % 

Max. 
Error % 

Mean 
Error % 

1A 

0.0010 

0.0005 

0.0010 

0.0005 

IB 

0.0100 

0.050 

0.0010 

0.0005 

1C 

1.3000 

0.6250 

0.0010 

0.0005 

2A 

0.0010 

0.0005 

0.0010 

0.0005 

2B 

0.0300 

0.0150 

0.0200 • 

0.0100 

2C 

0.3250 

0.1600 

0.2500 

0.1250 

3A 

0.0010 

0.0005 

0.0010 

0.0005 

3B 

0.0250 

0.0125 

0.0400 

0.0150 

3C 

0.2500 

0.1050 

0.1900 



0.0800 


Table 5.2 

Similarly as is done in the pure shear problem, a study of the 

f 

"assumed zero-stresses" is made for the pure bending problem. From the 
analytical solution, the stresses a v and T are zero. In the finite' 
element results, these stresses are not exactly zero, but are related 
to the element ratio and integration scheme used, among others. Fig- 
ures 5.18 to 5.23 present the variation of the shear stress x yv/ * along 
the y-axis, for the 3 maps and the two integration schemes studied. The 

choice of Tvu instead of a v is arbitrary since the variation pattern of 
xy x 

3 

both is essentially the same. An amplification factor of 10 is again 
used. It is observed that while the error in the 2x2 G.I, scheme 
monotonically decreases along the y-direction, for the 3x3 G.I. scheme, 
the error decreases harmonically in the same direction. Table 5.3 
presents a comparison of the "assumed-zero stress" error between the two 
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integration schemes studied. From the figures, it is observed that for 
cases A and B or any map, the difference between the two Integration 
schemes is minimal. However, for case C, specially Haps 2 and 3, 
there is a considerable improvement of the 3x3 G.I. over the 2x2 G.I. 
scheme. It Is interesting to note that the upper peaks of the curves 
of cases 2B, 2C, 3B, and 3C of the 3x3 G.I. scheme correspond to data 
provided by the same Gauss integration point of each element in the map. 



* 

The mean absolute error is being considered; that Is 

n 

% ME = l/n l X. 100 
i=l 1 x 

Table 5.3 

5.4 Execution Times 

As a final comparison between the tv/o integration schemes, the 
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execution times for each map obtained with the UNIVAC 1108 are studied. 

Since the only difference among the throe cases of any map is the 
relative size of the elements, an average time of the three cases of 
each map is taken. Figures 5.24 to 5.26 present the variation of 
execution time with map size for the two integration schemes studied 
and the three problems chosen. It is observed that for Maps 1 and 2 of 
the uniaxial tension and pure shear problems, the execution time is 
almost the same, but for Map 3 the difference is considerably larger. For 
the pure bending problem the pattern is different. Although for Map 1 
the time is the same, for Map 2, the time for the 3x3 G.I. scheme in- 
creases in 30%, but for Map 3 the difference is only of 17%. Note that 
the major difference in the three problems in terms of the finite element 
setup, is the prescription of the boundary conditions. Of course they 
affect the number of zero terms in the stiffness matrix and consequently the 
time of the routine that solves for the nodal displacements. From these 
time results we conclude that larger maps are necessary for better time 
comparisons between the two integration schemes. In [4] a time study of 
these integration procedures is presented for the cantilever beam problem. 
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6. CONCLUSIONS AND RECOMMENDATIONS 


It is concluded in the report that the eight-node isoparametric 
quadrilateral presents many of the features of the high-order elements 
with the advantage that there is no need to resort to expensive and 
complex integration schemes. 

It is shown that although the element contains in its formulation 
the exact solution for tension, shear and pure bending problems, the 
actual results are highly dependent upon the slenderness of the element 
in a given map, as well as the numerical Integration scheme used. As 
more terms of the displacement functions of the element enter in the solu- 
tion of a problem (linear terms for tension; quadratic terms for shear; 
cubic terms for pure bending, etc.) these effects must be carefully con- 
sidered in the accuracy of the solution. 

The results studied show that although' the 2x2 G.I. scheme provides 
adequate solutions, the penalty of using 3x3 G.I. Instead is not high 
in terms of storage and computation time, but the benefit is substantially 
higher. This' is particularly true in problems containing more complex 
geometry and loading, where the displacement solution is not contained in 
the element formulation. The same is valid for elasto-plastic problems 
where the location and size of the yielding zone is important, and where 
the construction of the stiffness matrix of an element that has yielded 
requires a higher-order integration procedure. 

In the transition to solve plasticity problems, apart from reserving 
some storage area for quantities such as accumulated forces and displace- 


ments, plastic strain energy, ‘plastic stresses and strains, the major 
change lies on the constitutive matrix, [D], This matrix is direct 
function of the state of loading, once yielding has occurred. This 
analysis is based on the theory of elasto-plastic flow described in 
detail in [8,9]. At any particular Gauss point where yielding has 
taken place the stress-strain matrix [D] must be updated at each loading 
step according to the current level of stress at that point. Hence, 
essentially only the assembly of the element: stiffness matrix gets 
changed since a different constitutive matrix might be needed for each 
Gauss point. 

Looking toward a general stress program for large computers, which 
would apply to realistic nonlinear material behavior of cracked structures, 
it is considered convenient to make the program capable to accept other 
element configurations, such as special crack-tip elements, or even 
other isoparametric elements. Furthermore, it should be possible to 
think in a quadrilateral having a cubic response along one side which 
contains four nodes, linear response along other side with only two corner 
nodes, and quadratic along the other two. The coupling of a higher order 
element to lower elements may be accomplished by constraining the coupling 
surface of the higher order element to displace in accordance to the lower 
order element. 

Finally, a desirable feature of any large computer program 'is an 
internal map generation algorithm. Because of the midside node in the 
isoparametric quadrilateral, the data preparation is considerably reduced 
by writing into the program an algorithm which interpolates the positions 
of midside node coordinates if the sides are straight. Only when the 


particular side is required to follow a curved boundary is it necessary 
to specify all intermediate nodes. However, algorithms for circum- 
ferential element profiles are easy to implement. 


ACKNOWLEDGMENTS 


The author is grateful to Professor J. L. Swedlow and Dr. J. R. Osi 
for their continuous advice and critical review of this paper. The 
assistance of Ms. Stella DeVito in the preparation of the manuscript is 
greatly appreciated. 


REFERENCES 


[13 Zienkiewicz, 0, C., The Finite Element Method in Engineering Science , 
McGraw-Hill , 1973. 

[2] Zienkiewicz, 0, C. t and C. J. Parekh, "Transient Field Problems: 

Two Dimensional and Three Dimensional Analysis by Isoparametric 
Finite Elements," IJNME , Vol. 2, 1970, pp. 61-71. 

[3] Ergatoudis, I., B. M. Irons, and 0. C. Zienkiewicz, "Curved, Iso- 
parametric Quadrilateral Elements for Finite Element Analysis," 
International journal of Solids Structures , Vol. 4, 1968, pp. 31-42. 

[4] Marino, C., "A Solution to the Cantilever Beam Problem using the 
Eight-Node Isoparametric Quadrilateral Report SM-75-3, Department 
of Mechanical Engineering, Carnegie-Mellon University, 1975. 

[5] Gupta, A. K., "Elasto-Plastic Analysis of Three-Dimensional Structures 
using the Isoparametric Element,'^ Nuclear Engineering and Design , 

Vol. 22, 1972, pp. 305-317, "" " ~ 

[6] Irons, B. M., "Engineering Applications of Numerical Integration in 
Stiffness Methods," A1AA Journal , Vol. 4, 1966, pp. 2035-2037. 

[7] Carnahan, B., H. A. Luther, and J, 0. Wilkes, Applied Numerical 
Methods , John Wiley & Sons, Inc., New York, 1969. 

[8] Swedlow, J. L., "A Procedure for Solving Problems of Elasto-Plastic 
Flow," Computer Structures , Vol. 3, 1973, pp. 879-898. 

[9] Swedlow, J. L., "Character of the Equations of Elasto-Plastic Flow 
in Three Independent Variables," International Journal of Non-Linear 
Mechanics , Vol. 3, 1968, pp. 325-336; Vol . 4, 1969, p. 77. 


APPENDIX A 


J 


The partial derivatives of the displacement functions for the 
8-Node Isoparametric quadrilateral are the following: 

(a) Corner Nodes 1 - 1,3, 5, 7 
In 5 : 

3Nj 

■gjT* = 1/4 5^0 + nn.j)(55.j + nn^ - 1) + 1/4(1 + £Cj)(l + nn^)c^ 

3N. 

Or, - 1/4 5^(1 + nn.j)(255. + nri^) (Al) 


In n: 


3N. 

= 1/4 (1 + 55^(5^ + nn i - 1) + 1/4(1 + 55.)(1 + nn i )n i 


8N. 

°r, = 1/4 n .(l + €e i ) Cee i + 2nn 1 ) (A2) 


(b) Midside Node, 5 ^ =0; i = 2,6 
In s: 



3N- 

g r = - 5(1 + nn i ) 

(A3) 

In n: 

3N • ,0 

ST= 1/2 n f O - £ 2 ) 

(A4) 

(c) Midside Nodes, 

J3 

II 

0 

— u 
II 

-Pi 

00 


In 5 : 

3N. ? 

it 1 " -!/ 2 M 1 

(A5) 

In n: 

3N. 

3n ° ' n(1 + “i 1 

(A 6 ) 


